module seq_domain_mct

  use shr_kind_mod, only: R8=>shr_kind_r8, IN=>shr_kind_in
  use shr_kind_mod, only: CL=>shr_kind_cl
  use shr_sys_mod,  only: shr_sys_flush, shr_sys_abort
  use shr_mpi_mod,  only: shr_mpi_min, shr_mpi_max

  use mct_mod
  use seq_cdata_mod
  use seq_comm_mct
  use seq_infodata_mod

  use map_atmlnd_mct
  use map_atmocn_mct
  use map_atmice_mct
  use map_iceocn_mct
  use map_snoglc_mct

  implicit none
  save
  private ! except

!--------------------------------------------------------------------------
! Public interfaces
!--------------------------------------------------------------------------

  public :: seq_domain_check_mct
  public :: domain_areafactinit_mct

!--------------------------------------------------------------------------
! Public variables
!--------------------------------------------------------------------------

  real(R8), parameter :: eps_tiny   = 1.0e-16_R8 ! roundoff eps
  real(R8), parameter :: eps_big    = 1.0e+02_R8 ! big eps
  real(R8), parameter :: eps_frac_samegrid = 1.0e-14_R8 ! epsilon for fractions for samegrid

!--------------------------------------------------------------------------
! Private interfaces
!--------------------------------------------------------------------------

#ifdef CPP_VECTOR
  logical :: usevector = .true.
#else
  logical :: usevector = .false.
#endif
  
#ifdef SYSUNICOS
  logical :: usealltoall = .true.
#else
  logical :: usealltoall = .false.
#endif
  
  private :: seq_domain_check_grid_mct

#include <mpif.h>  

!================================================================================
contains
!================================================================================

!================================================================================

  subroutine seq_domain_check_mct( cdata_a, cdata_i, cdata_l, cdata_o, cdata_r, cdata_g, cdata_s)

    !-----------------------------------------------------------
    !
    ! Arguments
    !
    type(seq_cdata), intent(in) :: cdata_a
    type(seq_cdata), intent(in) :: cdata_i
    type(seq_cdata), intent(in) :: cdata_l
    type(seq_cdata), intent(in) :: cdata_o
    type(seq_cdata), intent(in) :: cdata_r
    type(seq_cdata), intent(in) :: cdata_g
    type(seq_cdata), intent(in) :: cdata_s
    !
    ! Local variables
    !
    type(mct_gGrid)  , pointer :: atmdom_a   ! atm domain
    type(mct_gGrid)  , pointer :: icedom_i   ! ice domain
    type(mct_gGrid)  , pointer :: lnddom_l   ! lnd domain
    type(mct_gGrid)  , pointer :: ocndom_o   ! ocn domain
    type(mct_gGrid)  , pointer :: glcdom_g   ! glc domain
    type(mct_gGrid)  , pointer :: snodom_s   ! sno domain
    !
    type(mct_gsMap)  , pointer :: gsMap_a    ! atm global seg map 
    type(mct_gsMap)  , pointer :: gsMap_i    ! ice global seg map 
    type(mct_gsMap)  , pointer :: gsMap_l    ! lnd global seg map 
    type(mct_gsMap)  , pointer :: gsMap_o    ! ocn global seg map 
    type(mct_gsMap)  , pointer :: gsMap_r    ! ocn global seg map 
    type(mct_gsMap)  , pointer :: gsMap_g    ! glc global seg map 
    type(mct_gsMap)  , pointer :: gsMap_s    ! sno global seg map 
    !
    type(mct_gGrid) :: lnddom_a              ! lnd domain info on atm decomp
    type(mct_gGrid) :: icedom_a              ! ice domain info on atm decomp (all grids same)
    type(mct_gGrid) :: ocndom_a              ! ocn domain info on atm decomp (all grids same)
    type(mct_gGrid) :: icedom_o              ! ocn domain info on ocn decomp (atm/ocn grid different)
    type(mct_gGrid) :: snodom_g              ! sno domain info on glc decomp
    !
    integer(IN) :: mpicom_a                  ! atm mpicom
    integer(IN) :: mpicom_i                  ! ice mpicom
    integer(IN) :: mpicom_l                  ! lnd mpicom
    integer(IN) :: mpicom_o                  ! ocn mpicom
    integer(IN) :: mpicom_g                  ! glc mpicom
    integer(IN) :: mpicom_s                  ! sno mpicom
    !
    type(seq_infodata_type), pointer :: infodata
    !
    real(R8), pointer :: fracl(:)            ! land fraction on atm decomp 
    real(R8), pointer :: fraco(:)            ! ocn  fraction on atm decomp 
    real(R8), pointer :: fraci(:)            ! ice  fraction on atm decomp 
    real(R8), pointer :: maskl(:)            ! land mask on atm decomp (all grids same)
    real(R8), pointer :: maski(:)            ! ice  mask on atm decomp (all grids same)
    real(R8), pointer :: masko(:)            ! ocn  mask on atm decomp (all grids same)
    !
    integer(IN) :: n, kl, ko, ki             ! indicies
    integer(IN) :: k1,k2,k3                  ! indicies
    !
    logical      :: lnd_present              ! lnd present flag
    logical      :: ocn_present              ! ocn present flag
    logical      :: ice_present              ! ice present flag
    logical      :: glc_present              ! glc present flag
    logical      :: sno_present              ! sno present flag
    logical      :: rof_present              ! rof present flag
    logical      :: ocnrof_prognostic        ! ocn rof prognostic flag
    logical      :: samegrid_ao              ! atm ocn grid same
    logical      :: samegrid_ro              ! rof ocn grid same
    logical      :: samegrid_al              ! atm lnd grid same
    integer(IN)  :: rcode                    ! error status
    integer(IN)  :: atmsize                  ! local  size of atm  grid
    integer(IN)  :: lndsize                  ! local  size of land grid
    integer(IN)  :: ocnsize                  ! local  size of ocn  grid
    integer(IN)  :: icesize                  ! local  size of ice  grid
    integer(IN)  :: glcsize                  ! local  size of glc  grid
    integer(IN)  :: snosize                  ! local  size of sno  grid
    integer(IN)  :: gatmsize                 ! global size of atm  grid
    integer(IN)  :: glndsize                 ! global size of land grid
    integer(IN)  :: gocnsize                 ! global size of ocn  grid
    integer(IN)  :: grofsize                 ! global size of ocn  grid
    integer(IN)  :: gicesize                 ! global size of ice  grid
    integer(IN)  :: gglcsize                 ! global size of glc  grid
    integer(IN)  :: gsnosize                 ! global size of sno  grid
    integer(IN)  :: npts                     ! local size temporary
    integer(IN)  :: ier                      ! error code
    real(R8)     :: diff,dmaxo,dmaxi         ! difference tracker
    logical      :: iamroot                  ! local masterproc
    real(R8)     :: eps_frac                 ! epsilon for fractions
    real(R8)     :: eps_axmask               ! epsilon for masks, atm/lnd
    real(R8)     :: eps_axgrid               ! epsilon for grid coords, atm/lnd
    real(R8)     :: eps_axarea               ! epsilon for areas, atm/lnd
    real(R8)     :: eps_oimask               ! epsilon for masks, ocn/ice
    real(R8)     :: eps_oigrid               ! epsilon for grid coords, ocn/ice
    real(R8)     :: eps_oiarea               ! epsilon for areas, ocn/ice
    real(R8)     :: my_eps_frac              ! local eps_frac value
    real(r8)     :: rmin1,rmax1,rmin,rmax    ! local min max computation
    !
    real(R8),allocatable :: mask (:)         ! temporary real vector, domain mask
    !
    character(*),parameter :: F00 = "('(domain_check_mct) ',4a)"
    character(*),parameter :: F01 = "('(domain_check_mct) ',a,i6,a)"
    character(*),parameter :: F02 = "('(domain_check_mct) ',a,g23.15)"
    character(*),parameter :: F0R = "('(domain_check_mct) ',2A,2g23.15,A )"
    character(*),parameter :: subName = '(domain_check_mct) '
    !-----------------------------------------------------------

    call seq_comm_setptrs(CPLID,iamroot=iamroot)

    call seq_cdata_setptrs(cdata_a, infodata=infodata)
    call seq_infodata_GetData( infodata, samegrid_ao=samegrid_ao, samegrid_al=samegrid_al, &
         samegrid_ro=samegrid_ro)
    call seq_infodata_GetData( infodata, lnd_present=lnd_present, ocn_present=ocn_present, &
         ice_present=ice_present, glc_present=glc_present, sno_present=sno_present, &
         rof_present=rof_present, ocnrof_prognostic=ocnrof_prognostic)
    call seq_infodata_GetData( infodata, eps_frac=eps_frac, &
         eps_amask=eps_axmask, eps_agrid=eps_axgrid, eps_aarea=eps_axarea, &
         eps_omask=eps_oimask, eps_ogrid=eps_oigrid, eps_oarea=eps_oiarea )

    ! Get info

    call seq_cdata_setptrs(cdata_a, gsMap=gsMap_a, dom=atmdom_a, mpicom=mpicom_a)
    atmsize = mct_avect_lsize(atmdom_a%data)
    gatmsize = mct_gsMap_gsize(gsMap_a)

    if (lnd_present) then
       call seq_cdata_setptrs(cdata_l, gsMap=gsMap_l, dom=lnddom_l, mpicom=mpicom_l)
       lndsize = mct_avect_lsize(lnddom_l%data)
       glndsize = mct_gsMap_gsize(gsMap_l) 
       if (samegrid_al .and. gatmsize /= glndsize) then
          write(logunit,*) subname,' error: global atmsize = ',gatmsize,' global lndsize= ',glndsize
          call shr_sys_flush(logunit)
          call mct_die(subname,' atm and lnd grid must have the same global size')
       end if
       call mct_gGrid_init(oGGrid=lnddom_a, iGGrid=lnddom_l, lsize=atmsize)
       call mct_aVect_zero(lnddom_a%data)
       call map_lnd2atm_mct(cdata_l, lnddom_l%data, cdata_a, lnddom_a%data)
       allocate(maskl(atmsize),stat=rcode)
       if(rcode /= 0) call mct_die(subName,'allocate maskl')
       allocate(fracl(atmsize),stat=rcode)
       if(rcode /= 0) call mct_die(subName,'allocate fracl')
       call mct_aVect_exportRAttr(lnddom_a%data, 'mask', maskl, atmsize)
       call mct_aVect_exportRAttr(lnddom_a%data, 'frac', fracl, atmsize)
    endif

    if (ocn_present) then
       call seq_cdata_setptrs(cdata_o, gsMap=gsMap_o, dom=ocndom_o, mpicom=mpicom_o)
       ocnsize = mct_avect_lsize(ocndom_o%data)
       gocnsize = mct_gsMap_gsize(gsMap_o)
       if (samegrid_ao .and. gatmsize /= gocnsize) then
          write(logunit,*) subname,' error: global atmsize = ',gatmsize,' global ocnsize= ',gocnsize
          call shr_sys_flush(logunit)
          call mct_die(subname,' atm and ocn grid must have the same global size')
       end if
       call mct_gGrid_init(oGGrid=ocndom_a, iGGrid=ocndom_o, lsize=atmsize)
       call mct_aVect_zero(ocndom_a%data)
       call map_ocn2atm_mct(cdata_o, ocndom_o%data, cdata_a, ocndom_a%data)
       allocate(masko(atmsize),stat=rcode)
       if(rcode /= 0) call mct_die(subName,'allocate masko')
       allocate(fraco(atmsize),stat=rcode)
       if(rcode /= 0) call mct_die(subName,'allocate fraco')
       call mct_aVect_exportRAttr(ocndom_a%data, 'mask', masko, atmsize)
       if (samegrid_ao) then
          call mct_aVect_exportRattr(ocndom_a%data, 'frac', fraco, atmsize)
       else
          call mct_aVect_exportRattr(ocndom_a%data, 'mask', fraco, atmsize)
       endif
    endif
   
    if (ice_present) then
       call seq_cdata_setptrs(cdata_i, gsMap=gsMap_i, dom=icedom_i, mpicom=mpicom_i)
       icesize = mct_avect_lsize(icedom_i%data)
       gicesize = mct_gsMap_gsize(gsMap_i) 
       if (samegrid_ao .and. gatmsize /= gicesize) then
          write(logunit,*) subname,' error: global atmsize = ',gatmsize,' global icesize= ',gicesize
          call shr_sys_flush(logunit)
          call mct_die(subname,' atm and ice grid must have the same global size')
       end if
       call mct_gGrid_init(oGGrid=icedom_a, iGGrid=icedom_i, lsize=atmsize)
       call mct_aVect_zero(icedom_a%data)
       call map_ice2atm_mct(cdata_i, icedom_i%data, cdata_a, icedom_a%data)
       allocate(maski(atmsize),stat=rcode)
       if(rcode /= 0) call mct_die(subName,'allocate maski')
       allocate(fraci(atmsize),stat=rcode)
       if(rcode /= 0) call mct_die(subName,'allocate fraci')
       call mct_aVect_exportRAttr(icedom_a%data, 'mask', maski, atmsize)
       if (samegrid_ao) then
          call mct_aVect_exportRattr(icedom_a%data, 'frac', fraci, atmsize)
       else
          call mct_aVect_exportRattr(icedom_a%data, 'mask', fraci, atmsize)
       endif
    endif

    if (sno_present .and. glc_present) then
       call seq_cdata_setptrs(cdata_g, gsMap=gsMap_g, dom=glcdom_g, mpicom=mpicom_g)
       glcsize = mct_avect_lsize(glcdom_g%data)
       gglcsize = mct_gsMap_gsize(gsMap_g) 
       call seq_cdata_setptrs(cdata_s, gsMap=gsMap_s, dom=snodom_s, mpicom=mpicom_s)
       snosize = mct_avect_lsize(snodom_s%data)
       gsnosize = mct_gsMap_gsize(gsMap_s) 
       if (gglcsize /= gsnosize) then
          write(logunit,*) subname,' error: global glcsize = ',gglcsize,' global snosize= ',gsnosize
          call shr_sys_flush(logunit)
          call mct_die(subname,' glc and sno grid must have the same global size')
       end if
       call mct_gGrid_init(oGGrid=snodom_g, iGGrid=snodom_s, lsize=glcsize)
       call mct_aVect_zero(snodom_g%data)
       call map_sno2glc_mct(cdata_s, snodom_s%data, cdata_g, snodom_g%data)
       if (iamroot) write(logunit,F00) ' --- checking glc/sno domains ---'
       npts = glcsize
       allocate(mask(npts),stat=rcode)
       if(rcode /= 0) call mct_die(subName,'allocate mask')
       call mct_aVect_getRAttr(snodom_g%data,"mask",mask,rcode)
       where (mask < eps_axmask) mask = 0.0_r8
       call seq_domain_check_grid_mct(glcdom_g%data, snodom_g%data, 'mask', eps=eps_axmask, mpicom=mpicom_g, mask=mask)
       call seq_domain_check_grid_mct(glcdom_g%data, snodom_g%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_g, mask=mask)
       call seq_domain_check_grid_mct(glcdom_g%data, snodom_g%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_g, mask=mask)
       call seq_domain_check_grid_mct(glcdom_g%data, snodom_g%data, 'area', eps=eps_axarea, mpicom=mpicom_g, mask=mask)
       deallocate(mask,stat=rcode)
       if(rcode /= 0) call mct_die(subName,'deallocate mask')
    endif

    if (ice_present .and. ocn_present) then
       if (gocnsize /= gicesize) then
          write(logunit,*) subname,' error: global ocnsize = ',gocnsize,' global icesize= ',gicesize
          call shr_sys_flush(logunit)
          call mct_die(subname,' ocean and ice grid must have the same global size')
       endif
       call mct_gGrid_init(oGGrid=icedom_o, iGGrid=icedom_i, lsize=ocnsize)
       call mct_aVect_zero(icedom_o%data)
       call map_ice2ocn_mct(cdata_i, icedom_i%data, cdata_o, icedom_o%data)
    end if

    if (rof_present .and. ocnrof_prognostic .and. samegrid_ro) then
       call seq_cdata_setptrs(cdata_r, gsMap=gsMap_r)
       grofsize = mct_gsMap_gsize(gsMap_r)
       if (gocnsize /= grofsize) then
          write(logunit,*) subname,' error: global ocnsize = ',gocnsize,' global rofsize= ',grofsize
          call shr_sys_flush(logunit)
          call mct_die(subname,' ocean and rof grid must have the same global size')
       endif
    end if

    !------------------------------------------------------------------------------
    ! Check ice/ocean grid consistency
    !------------------------------------------------------------------------------

     if (ocn_present .and. ice_present) then
!    if (samegrid_oi) then       ! doesn't yet exist

       npts = ocnsize
       allocate(mask(npts),stat=rcode)
       if(rcode /= 0) call mct_die(subName,'allocate mask')

       if (iamroot) write(logunit,F00) ' --- checking ocn/ice domains ---'
       call seq_domain_check_grid_mct(ocndom_o%data, icedom_o%data,'mask', eps=eps_oigrid, mpicom=mpicom_o)
       call mct_aVect_getRAttr(ocndom_o%data,"mask",mask,rcode)
       where (mask < eps_oimask) mask = 0.0_r8

       call seq_domain_check_grid_mct(ocndom_o%data, icedom_o%data,'lat' , eps=eps_oigrid, mpicom=mpicom_o, mask=mask)
       call seq_domain_check_grid_mct(ocndom_o%data, icedom_o%data,'lon' , eps=eps_oigrid, mpicom=mpicom_o, mask=mask)
       call seq_domain_check_grid_mct(ocndom_o%data, icedom_o%data,'area', eps=eps_oiarea, mpicom=mpicom_o, mask=mask)

       deallocate(mask,stat=rcode)
       if(rcode /= 0) call mct_die(subName,'deallocate mask')

!    endif
     endif

    !------------------------------------------------------------------------------
    ! Check atm/lnd grid consistency
    !------------------------------------------------------------------------------

    if (lnd_present .and. samegrid_al) then
       if (iamroot) write(logunit,F00) ' --- checking atm/land domains ---'
       call seq_domain_check_grid_mct(atmdom_a%data, lnddom_a%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_a, mask=maskl)
       call seq_domain_check_grid_mct(atmdom_a%data, lnddom_a%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_a, mask=maskl)
       call seq_domain_check_grid_mct(atmdom_a%data, lnddom_a%data, 'area', eps=eps_axarea, mpicom=mpicom_a, mask=maskl)
    endif

    !------------------------------------------------------------------------------
    ! Check atm/ocn and atm/ice grid consistency (if samegrid)
    !------------------------------------------------------------------------------

    if (ice_present .and. samegrid_ao) then
       if (iamroot) write(logunit,F00) ' --- checking atm/ice domains ---'
       call seq_domain_check_grid_mct(atmdom_a%data, icedom_a%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_a, mask=maski)
       call seq_domain_check_grid_mct(atmdom_a%data, icedom_a%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_a, mask=maski)
       call seq_domain_check_grid_mct(atmdom_a%data, icedom_a%data, 'area', eps=eps_axarea, mpicom=mpicom_a, mask=maski)
    endif

    if (ocn_present .and. samegrid_ao) then
       if (iamroot) write(logunit,F00) ' --- checking atm/ocn domains ---'
       call seq_domain_check_grid_mct(atmdom_a%data, ocndom_a%data, 'lat' , eps=eps_axgrid, mpicom=mpicom_a, mask=masko)
       call seq_domain_check_grid_mct(atmdom_a%data, ocndom_a%data, 'lon' , eps=eps_axgrid, mpicom=mpicom_a, mask=masko)
       call seq_domain_check_grid_mct(atmdom_a%data, ocndom_a%data, 'area', eps=eps_axarea, mpicom=mpicom_a, mask=masko)
    endif

    !------------------------------------------------------------------------------
    ! Check consistency of land fraction with ocean mask on grid
    !------------------------------------------------------------------------------

    my_eps_frac = eps_frac
    if (samegrid_ao) my_eps_frac = eps_frac_samegrid
    if (.not. samegrid_al) my_eps_frac = eps_big

    if (iamroot) write(logunit,F00) ' --- checking fractions in domains ---'
    dmaxi = 0.0_R8
    dmaxo = 0.0_R8
    do n = 1,atmsize
       if (lnd_present .and. ice_present) then
          diff = abs(1._r8 - fracl(n) - fraci(n))
          dmaxi = max(diff,dmaxi)
          if (diff > my_eps_frac) then
             write(logunit,*)'inconsistency between land fraction and sea ice fraction'
             write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraci= ',fraci(n),' sum= ',fracl(n)+fraci(n)
             call shr_sys_flush(logunit)
             call mct_die(subname,'inconsistency between land fraction and sea ice fraction')
          end if
          if ((1._r8-fraci(n)) > eps_frac .and. fracl(n) < eps_tiny) then
             write(logunit,*)'inconsistency between land mask and sea ice mask'
             write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraci= ',fraci(n)
             call shr_sys_flush(logunit)
             call mct_die(subname,' inconsistency between land mask and sea ice mask')
          end if
       endif
       if (lnd_present .and. ocn_present) then
          diff = abs(1._r8 - fracl(n) - fraco(n))
          dmaxo = max(diff,dmaxo)
          if (diff > my_eps_frac) then
             write(logunit,*)'inconsistency between land fraction and ocn land fraction'
             write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraco= ',fraco(n),' sum= ',fracl(n)+fraco(n)
             call shr_sys_flush(logunit)
             call mct_die(subname,' inconsistency between land fraction and ocn land fraction')
          end if
          if ((1._r8-fraco(n)) > eps_frac .and. fracl(n) < eps_tiny) then
             write(logunit,*)'inconsistency between land mask and ocn land mask'
             write(logunit,*)'n= ',n,' fracl= ',fracl(n),' fraco= ',fraco(n)
             call shr_sys_flush(logunit)
             call mct_die(subname,' inconsistency between land mask and ocn land mask')
          end if
       endif
    end do 
    if (iamroot) then
       write(logunit,F02) ' maximum           difference for ofrac sum ',dmaxo
       write(logunit,F02) ' maximum           difference for ifrac sum ',dmaxi
       write(logunit,F02) ' maximum allowable difference for  frac sum ',my_eps_frac
       write(logunit,F02) ' maximum allowable tolerance for valid frac ',eps_frac
       call shr_sys_flush(logunit)
    endif

    !------------------------------------------------------------------------------
    ! Set atm and lnd ascale
    !------------------------------------------------------------------------------
    if (lnd_present .and. ocn_present) then
       k1 = mct_aVect_indexRa(atmdom_a%data,"ascale",perrWith='domain_check ascale')
       do k2 = 1,atmsize
          if (fracl(k2) /= 0.0_r8) then
             atmdom_a%data%rAttr(k1,k2) = (1.0_r8-fraco(k2))/(fracl(k2))
          else
             atmdom_a%data%rAttr(k1,k2) = 0.0
          endif
       enddo
       call map_atm2lnd_mct(cdata_a, atmdom_a%data, cdata_l, lnddom_l%data, fluxlist='ascale')
    elseif (lnd_present .and. ice_present) then
       k1 = mct_aVect_indexRa(atmdom_a%data,"ascale",perrWith='domain_check ascale')
       do k2 = 1,atmsize
          if (fracl(k2) /= 0.0_r8) then
             atmdom_a%data%rAttr(k1,k2) = (1.0_r8-fraci(k2))/(fracl(k2))
          else
             atmdom_a%data%rAttr(k1,k2) = 0.0
          endif
       enddo
       call map_atm2lnd_mct(cdata_a, atmdom_a%data, cdata_l, lnddom_l%data, fluxlist='ascale')
    endif

    if (lnd_present .and. (ocn_present .or. ice_present)) then
       k1 = mct_aVect_indexRa(atmdom_a%data,"ascale",perrWith='domain_check atm ascale')
       rmin1 = minval(atmdom_a%data%rAttr(k1,:))
       rmax1 = maxval(atmdom_a%data%rAttr(k1,:))
       call shr_mpi_min(rmin1,rmin,mpicom_a)
       call shr_mpi_max(rmax1,rmax,mpicom_a)
       if (iamroot) write(logunit,F0R) trim(subname),' : min/max ascale ',rmin,rmax,' atmdom_a'

       k1 = mct_aVect_indexRa(lnddom_l%data,"ascale",perrWith='domain_check lnd ascale')
       rmin1 = minval(lnddom_l%data%rAttr(k1,:))
       rmax1 = maxval(lnddom_l%data%rAttr(k1,:))
       call shr_mpi_min(rmin1,rmin,mpicom_l)
       call shr_mpi_max(rmax1,rmax,mpicom_l)
       if (iamroot) write(logunit,F0R) trim(subname),' : min/max ascale ',rmin,rmax,' lnddom_l'
    endif

    !------------------------------------------------------------------------------
    ! Clean up allocated memory
    !------------------------------------------------------------------------------

    if (lnd_present) then
       deallocate(fracl,stat=rcode)
       if(rcode /= 0) call mct_die(subName,'deallocate fracl')
       deallocate(maskl,stat=rcode)
       if(rcode /= 0) call mct_die(subName,'deallocate maskl')
       call mct_gGrid_clean(lnddom_a, rcode)
       if(rcode /= 0) call mct_die(subName,'clean lnddom_a')
    endif

    if (ocn_present) then
       deallocate(fraco,stat=rcode)
       if(rcode /= 0) call mct_die(subName,'deallocate fraco')
       deallocate(masko,stat=rcode)
       if(rcode /= 0) call mct_die(subName,'deallocate masko')
       call mct_gGrid_clean(ocndom_a, rcode)
       if(rcode /= 0) call mct_die(subName,'clean ocndom_a')
    endif

    if (ice_present) then
       deallocate(fraci,stat=rcode)
       if(rcode /= 0) call mct_die(subName,'deallocate fraci')
       deallocate(maski,stat=rcode)
       if(rcode /= 0) call mct_die(subName,'deallocate maski')
       call mct_gGrid_clean(icedom_a, rcode)
       if(rcode /= 0) call mct_die(subName,'clean icedom_o')
    endif

    if (ocn_present .and. ice_present) then
       call mct_gGrid_clean(icedom_o, rcode)
       if(rcode /= 0) call mct_die(subName,'clean icedom_o')
    endif

  end subroutine seq_domain_check_mct

!===============================================================================
  
  subroutine seq_domain_check_grid_mct(dom1, dom2, attr, eps, mpicom, mask)
   
    !-----------------------------------------------------------

    ! Arguments

    type(mct_aVect) , intent(in) :: dom1
    type(mct_aVect) , intent(in) :: dom2
    character(len=*), intent(in) :: attr   ! grid attribute to compare
    real(R8)        , intent(in) :: eps    ! error condition for compare
    integer(IN)     , intent(in) :: mpicom
    real(R8)        , intent(in), optional :: mask(:)

    ! Local variables

    integer(in)       :: n,ndiff            ! indices
    integer(in)       :: npts1,npts2,npts   ! counters
    integer(in)       :: rcode              ! error code
    real(R8)          :: diff,max_diff      ! temporaries
    real(R8)          :: tot_diff           ! maximum diff across all pes
    integer(IN)       :: ier                ! error code
    real(R8), pointer :: data1(:)           ! temporaries
    real(R8), pointer :: data2(:)           ! temporaries
    real(R8), pointer :: lmask(:)           ! temporaries
    logical           :: iamroot            ! local masterproc

    character(*),parameter :: F00 = "('(domain_check_grid_mct) ',4a)"
    character(*),parameter :: F01 = "('(domain_check_grid_mct) ',a,i12,a)"
    character(*),parameter :: F02 = "('(domain_check_grid_mct) ',2a,g23.15)"
    character(*),parameter :: F0R = "('(domain_check_grid_mct) ',2A,2g23.15,A )"
    character(*),parameter :: subName = '(domain_check_grid_mct) '
    !-----------------------------------------------------------

    call seq_comm_setptrs(CPLID,iamroot=iamroot)

   npts1 = mct_aVect_lsize(dom1)
   npts2 = mct_aVect_lsize(dom2)
   npts  = npts1
   
   if (npts1 == npts2) then
      if (iamroot) write(logunit,F01) " the domain size is = ", npts
   else
      write(logunit,*) trim(subname)," domain size #1 = ", npts1
      write(logunit,*) trim(subname)," domain size #2 = ", npts2
      write(logunit,*) trim(subname)," ERROR: domain size mis-match"
      call mct_die(subName,"ERROR: domain size mis-match")
   end if
   
   allocate(data1(npts),stat=rcode)
   if(rcode /= 0) call mct_die(subName,'allocate data1')
   allocate(data2(npts),stat=rcode)
   if(rcode /= 0) call mct_die(subName,'allocate data2')
   allocate(lmask(npts),stat=rcode)
   if(rcode /= 0) call mct_die(subName,'allocate lmask')

   call mct_aVect_exportRAttr(dom1, trim(attr), data1, npts)
   call mct_aVect_exportRAttr(dom2, trim(attr), data2, npts)
   lmask = 1.0_r8
   if (present(mask)) then
      if (size(mask) /= npts) then
         call mct_die(subName,"ERROR: mask size mis-match")
      endif
      lmask = mask
   endif

   ! --- adjust lons to address wraparound issues, we're assuming degree here! ---

   if (trim(attr) == "lon") then
      do n = 1,npts
         if (data2(n) > data1(n)) then
            do while ( (data1(n)+360.0_R8) < (data2(n)+180.0_R8) ) ! longitude is periodic
               data1(n) = data1(n) + 360.0_R8
            end do
         else
            do while ( (data2(n)+360.0_R8) < (data1(n)+180.0_R8) ) ! longitude is periodic
               data2(n) = data2(n) + 360.0_R8
            end do
         endif
      enddo
   endif

   ! Only check consistency where mask is greater than zero, if mask is present

   max_diff = 0.0_R8
   ndiff = 0
   do n=1,npts
      if (lmask(n) > eps_tiny) then
         diff = abs(data1(n)-data2(n))
         max_diff = max(max_diff,diff)
         if (diff > eps) then
!debug            write(logunit,*)'n= ',n,' data1= ',data1(n),' data2= ',data2(n),' diff= ',diff, ' eps= ',eps
            ndiff = ndiff + 1
         endif
      end if
   end do

   call mpi_reduce(max_diff,tot_diff,1,MPI_REAL8,MPI_MAX,0,mpicom,ier)
   if (iamroot) then
      write(logunit,F02) " maximum           difference for ",trim(attr),tot_diff
      write(logunit,F02) " maximum allowable difference for ",trim(attr),eps
      call shr_sys_flush(logunit)
   endif
   call mpi_barrier(mpicom,ier)

   if (ndiff > 0) then
      write(logunit,*) trim(subname)," ERROR: incompatible domain grid coordinates"
      call shr_sys_flush(logunit)
      call mct_die(subName,"incompatible domain grid coordinates")
   endif
   
   deallocate(data1,stat=rcode)
   if(rcode /= 0) call mct_die(subName,'deallocate data1')
   deallocate(data2,stat=rcode)
   if(rcode /= 0) call mct_die(subName,'deallocate data2')
   deallocate(lmask,stat=rcode)
   if(rcode /= 0) call mct_die(subName,'deallocate lmask')

 end subroutine seq_domain_check_grid_mct

!===============================================================================

 subroutine domain_areafactinit_mct( cdata, mdl2drv, drv2mdl, comment)
    !-----------------------------------------------------------
    !
    ! Arguments
    !
    type(seq_cdata) , intent(inout) :: cdata
    real(R8),pointer                :: mdl2drv(:)
    real(R8),pointer                :: drv2mdl(:)
    character(len=*),optional,intent(in) :: comment
    !
    ! Local variables
    !
    type(seq_infodata_type),pointer :: infodata
    type(mct_gGrid),pointer:: domain
    integer                :: ID
    integer                :: mpicom
    logical                :: iamroot
    logical                :: samegrid_ao, samegrid_al
    integer                :: j1,j2,m1,n,rcode
    integer                :: gridsize,m2dsize,d2msize
    real(r8)               :: rmin1,rmax1,rmin,rmax
    real(r8)               :: rmask,rarea,raream
    character(cl)          :: lcomment
    character(len=*),parameter :: subName = '(domain_areafactinit_mct) '
    character(len=*),parameter :: F0R = "(2A,2g23.15,A )"
    !
    !-----------------------------------------------------------

    lcomment = ''
    if (present(comment)) lcomment = comment

    call seq_cdata_setptrs(cdata, ID=ID, dom=domain, infodata=infodata)
    call seq_comm_setptrs(ID, mpicom=mpicom, iamroot=iamroot) 
    call seq_infodata_GetData( infodata, samegrid_ao=samegrid_ao, samegrid_al=samegrid_al)

    ! get sizes

    gridsize = mct_gGrid_lsize(domain)
    allocate(drv2mdl(gridsize),mdl2drv(gridsize),stat=rcode)
    if(rcode /= 0) call mct_die(subname,'allocate area correction factors')

    j1 = mct_gGrid_indexRA(domain,"area"    ,dieWith=subName)
    j2 = mct_gGrid_indexRA(domain,"aream"   ,dieWith=subName)
    m1 = mct_gGrid_indexRA(domain,"mask"    ,dieWith=subName)

    mdl2drv(:)=1.0_R8
    drv2mdl(:)=1.0_R8

    if (samegrid_ao .and. samegrid_al) then
        ! default 1.0
    else
       do n=1,gridsize
          rmask  = domain%data%rAttr(m1,n)
          rarea  = domain%data%rAttr(j1,n)
          raream = domain%data%rAttr(j2,n)
          ! lihuimin, TODO, FOR DEBUG
          !if ( abs(rmask) >= 1.0e-06) then
          !  if (rarea * raream /= 0.0_r8) then
          !     mdl2drv(n) = rarea/raream
          !     drv2mdl(n) = 1.0_R8/mdl2drv(n)
          !     !if (mdl2drv(n) > 10.0 .or. mdl2drv(n) < 0.1) then
          !     !   write(logunit,*) trim(subname),' WARNING area,aream= ', &
          !     !      domain%data%rAttr(j1,n),domain%data%rAttr(j2,n),' in ',n,gridsize
          !     !endif
          !  else
          !     write(logunit,*) trim(subname),' ERROR area,aream= ', &
          !        rarea,raream,' in ',n,gridsize
          !     call shr_sys_flush(logunit)
          !     call shr_sys_abort()
          !  endif
          !endif
       enddo
    end if
       
    rmin1 = minval(mdl2drv)
    rmax1 = maxval(mdl2drv)
    call shr_mpi_min(rmin1,rmin,mpicom)
    call shr_mpi_max(rmax1,rmax,mpicom)
    if (iamroot) write(logunit,F0R) trim(subname),' : min/max mdl2drv ',rmin,rmax,trim(lcomment)

    rmin1 = minval(drv2mdl)
    rmax1 = maxval(drv2mdl)
    call shr_mpi_min(rmin1,rmin,mpicom)
    call shr_mpi_max(rmax1,rmax,mpicom)
    if (iamroot) write(logunit,F0R) trim(subname),' : min/max drv2mdl ',rmin,rmax,trim(lcomment)

 end subroutine domain_areafactinit_mct

!===============================================================================

end module seq_domain_mct



